## Additional Attainment Status Analysis (Results are estimated by lfe package version 2.8-7)
## These results are reported in Appendix Table A8-A9
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(readxl)
library('usmap')
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
## Generate county info
data3$StateName=as.character(data3$StateName)
data3$County=as.character(data3$County)
county_list=unique(data3[c("StateName","County")])
## Fix county names to match non-attainment data
county_list$County1=county_list$County
county_list=subset(county_list, !(is.na(County1)))
county_list$County1[which(county_list$County1=="Miami Dade")]="Miami-Dade"
county_list$County1[which(county_list$County1=="St Lucie")]="St. Lucie"
county_list$County1[which(county_list$StateName=="Louisiana")]=paste(county_list$County1[which(county_list$StateName=="Louisiana")], ' Parish', sep='')
county_list$County1[which(county_list$County1=="St Mary Parish")]="St. Mary Parish"
county_list$County1[which(county_list$County1=="DeSoto Parish")]="De soto Parish"
county_list$County1[which(county_list$County1=="St James Parish")]="St. James Parish"
county_list$County1[which(county_list$County1=="St Charles Parish")]="St. Charles Parish"
county_list$County1[which(county_list$County1=="St Bernard Parish")]="St. Bernard Parish"
county_list$County1[which(county_list$County1=="Queen Annes")]=paste("Queen Anne's")
county_list$County1[which(county_list$County1=="Prince Georges")]=paste("Prince George's")
county_list$County1[which(county_list$County1=="St Clair")]="St. Clair"
county_list$County1[which(county_list$County1=="St Francois")]="St. Francois"
county_list$County1[which(county_list$County1=="St Louis")]="St. Louis"
county_list$County1[which(county_list$County1=="St Charles")]="St. Charles"
county_list$County1[which(county_list$County1=="St Louis City")]="St. Louis City"
county_list$County1[which(county_list$County1=="DeWitt")]="De Witt"
county_list$County1[which(county_list$County1=="Dewitt")]="De Witt"
county_list$County1[which(county_list$County1=="St Joseph")]="St. Joseph"
county_list$County1[which(county_list$County1=="St Croix")]="St. Croix"
county_list$County1[which(county_list$County1=="St Lawrence")]="St. Lawrence"
county_list$fips=NA
for (i in 1:nrow(county_list)) {
  county_list$fips[i]=fips(county_list$StateName[i], county_list$County1[i])
  print(i)
}
county_list=county_list[c("StateName","County","fips")]

data3=left_join(data3, county_list, by=c("StateName","County"))
data3$fips=as.numeric(data3$fips)

########################
## Extract attainment information
nonattain=read_xls("raw data/nonattain.xls") ## From EPA NAAQS nonattainment status
nonattain$fips=as.numeric(nonattain$fips_state)*1000+as.numeric(nonattain$fips_cnty)
nonattain=nonattain[c('fips',"pollutant","pw_2013","pw_2014","pw_2015","pw_2016","pw_2017","pw_2018")]
nonattain$Year_2013="Attain"
nonattain$Year_2013[which(nonattain$pw_2013=='P')]='Partial'
nonattain$Year_2013[which(nonattain$pw_2013=='W')]='Whole'
nonattain$Year_2014="Attain"
nonattain$Year_2014[which(nonattain$pw_2014=='P')]='Partial'
nonattain$Year_2014[which(nonattain$pw_2014=='W')]='Whole'
nonattain$Year_2015="Attain"
nonattain$Year_2015[which(nonattain$pw_2015=='P')]='Partial'
nonattain$Year_2015[which(nonattain$pw_2015=='W')]='Whole'
nonattain$Year_2016="Attain"
nonattain$Year_2016[which(nonattain$pw_2016=='P')]='Partial'
nonattain$Year_2016[which(nonattain$pw_2016=='W')]='Whole'
nonattain$Year_2017="Attain"
nonattain$Year_2017[which(nonattain$pw_2017=='P')]='Partial'
nonattain$Year_2017[which(nonattain$pw_2017=='W')]='Whole'
nonattain$Year_2018="Attain"
nonattain$Year_2018[which(nonattain$pw_2018=='P')]='Partial'
nonattain$Year_2018[which(nonattain$pw_2018=='W')]='Whole'


nonattain=subset(nonattain, Year_2018!="Attain"|Year_2017!="Attain"|Year_2016!="Attain"
                 |Year_2015!="Attain"|Year_2014!="Attain"|Year_2013!="Attain")

nonattain_all=subset(nonattain, Year_2018!="Attain"&Year_2017!="Attain"&Year_2016!="Attain"
                     &Year_2015!="Attain"&Year_2014!="Attain"&Year_2013!="Attain")

nonattain_2018=unique(nonattain$fips[which(nonattain$Year_2018!="Attain")])   

nonattain=unique(nonattain$fips)

nonattain_all=unique(nonattain_all$fips)

switcher=subset(nonattain, !(nonattain%in%nonattain_all))

switcher_non2018=subset(switcher, switcher%in%nonattain_2018)
switcher_attain2018=subset(switcher, !(switcher%in%nonattain_2018))

data3$attain=1
data3$attain[which(data3$fips%in%nonattain)]=0
data3$attain[which(data3$fips%in%switcher_non2018)]=2
data3$attain[which(data3$fips%in%switcher_attain2018)]=3

data3$switching_PM=0
data3$switching_PM[which(data3$fips%in%c(13015,13115,13207,18043,18077,21015,29071,
                                       29183,39025,42021,47001,47145,55079,54079))]=1


data3$shutdown_attain=0
data3$shutdown_attain[which(data3$shutdown==1&data3$attain==1)]=1

data3$shutdown_nonattain=0
data3$shutdown_nonattain[which(data3$shutdown==1&data3$attain==0)]=1

data4=data3
data4$shutdown_attain=0
data4$shutdown_attain[which(data4$shutdown==1&data4$attain==1)]=1

data4$shutdown_nonattain=0
data4$shutdown_nonattain[which(data4$shutdown==1&data4$attain==0)]=1

data4$shutdown_switch=0
data4$shutdown_switch[which(data4$shutdown==1&data4$attain>=2)]=1

data4$shutdown_switch_non2018=0
data4$shutdown_switch_non2018[which(data4$shutdown==1&data4$attain==2)]=1


data4$shutdown_switch_attain2018=0
data4$shutdown_switch_attain2018[which(data4$shutdown==1&data4$attain==3)]=1

data4$shutdown_switch_PM=0
data4$shutdown_switch_PM[which(data4$shutdown==1&data4$switching_PM==1)]=1

data4$shutdown_switch_nonPM=0
data4$shutdown_switch_nonPM[which(data4$shutdown==1&data4$switching_PM!=1)]=1

data5=subset(data4, T>40&T<127)




##############################################################################################################################
##############################################################################################################################
## Baseline Model
data_base1=subset(data4, T<99)
data_base2=subset(data5, T<99)
## Overall effect

all_AOD_long=felm(AOD~ shutdown_attain + shutdown_nonattain + shutdown_switch + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data_base1)
summary(all_AOD_long)

all_AOD_long2=felm(AOD~ shutdown_attain + shutdown_nonattain  + shutdown_switch_non2018 + shutdown_switch_attain2018  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data_base1)
summary(all_AOD_long2)

all_AOD_short=felm(AOD~ shutdown_attain + shutdown_nonattain + shutdown_switch  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data_base2)
summary(all_AOD_short)

all_AOD_short2=felm(AOD~ shutdown_attain + shutdown_nonattain + shutdown_switch_non2018 + shutdown_switch_attain2018  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                    data = data_base2)
summary(all_AOD_short2)


stargazer(all_AOD_long2,all_AOD_short2, digits=3)

stargazer(all_AOD_long,all_AOD_short, digits=3)


